1. Field of the Invention
The present invention relates to a simulation method for obtaining a physical quantity such as a potential inside an electrical element such as a semiconductor device, or parasitic element constants such as capacitance, resistance and inductance by numerical analysis.
2. Description of the Background Art
As an electrical element such as a semiconductor device becomes smaller, an effect that a parasitic element produces on the electrical element becomes relatively greater, so a quantitative grasp on the effect is required. It is generally difficult to obtain the parasitic element constant which quantitatively indicates the effect of the parasitic element by measurement of a practical device, and therefore the parasitic element constant is obtained by a simulation of numerical analysis with a computer.
Through this simulation, a numerical analysis is performed on a physical quantity (e.g., potential, temperature and the like) in a given position inside the electrical element, accompanying computation of the parasitic element constant. As a method of numerical analysis, the FEM (finite element method), the FDM (finite difference method) or the BEM (boundary element method) is often used.
Discussion will be made below on an exemplary case where a value of parasitic capacitance of a conductor in a two-dimensional analysis region is numerically analyzed by the finite difference method using the control volume method.
An direct object of the finite difference method is to obtain a physical quantity in a point of intersection (referred to as xe2x80x9cnodexe2x80x9d) of mesh (hereinafter, for simple discussion, a potential is adopted as a specific example of the physical quantity, and discussion is based thereon unless otherwise specified).
FIG. 9 illustrates a structure of a region to be analyzed. This figure shows a partial region VP of an electrical element and the region VP consists of a conductor region CR1 and a dielectric region DR1 having a dielectric constant of ∈. The region VP is divided into meshes on expression (the conductor region CR1 is not divided into meshes since uniform potential exists in an electrostatic field inside the conductor region CR1). The conductor region CR1 has material boundaries C1 and C2 and is separated thereby from others. The dielectric region DR1 has material boundaries C2 and C3 and analysis boundaries D1 an D2 and is separated thereby from others. In the material boundary C3, the dielectric region DR1 adjoins a not-shown conductor region. Herein, the analysis boundary refers to a limit of an analyzable region and the material boundary refers to a boundary where different materials adjoin.
It is assumed in FIG. 9 that for example, the potentials at the nodes on the material boundaries C2 and C3 and the normal components of electric fields (derivative of potential to position) relative to interfaces (hereinafter, referred to as xe2x80x9cnormal differential of potentialxe2x80x9d) at the nodes on the analysis boundaries D1 and D2 are given as boundary conditions. A node having an unknown quantity to be obtained is represented by a blank circle (◯), a node given a value of potential as a boundary condition is represented by a solid circle (xe2x97xaf) and a node given a value of normal differential of potential as a boundary condition is represented by a double circle ((⊚). In a case of FIG. 9, the unknown quantity to be obtained is a potential at the node represented by ◯ and there are 36 nodes of ◯ in the whole region VP. Values of potentials at the 36 points should be obtained under the above boundary conditions.
From the potential distribution thus obtained, the amount of electric charges that the conductor region CR1 and other not-shown conductor regions have under the respective boundary conditions is calculated, to obtain the parasitic capacitance of the conductor through numerical analysis.
As another exemplary case for obtaining the parasitic capacitance of the conductor in the two-dimensional analysis region, numeric analysis by the boundary element method will be discussed below.
FIG. 10 illustrates a structure of a region to be analyzed. This figure shows the same region VP as FIG. 9 by different expression. In the boundary element method, insince the boundary is expressed like a polygonal line and an object is to obtain the potential or its normal differential at each segment line (hereinafter, referred to as xe2x80x9cboundary elementxe2x80x9d), the inside of the dielectric region DR1 is not expressed in a form of mesh. In FIG. 10, the boundary elements are represented by the nodes (xe2x97xaf, ⊚) which are points of intersection of the boundary elements, like in FIG. 9.
In the boundary element method, at each boundary element, the value of either one of the potential and its normal differential is given as a boundary condition and the value of the other is the unknown quantity to be obtained. Therefore, in the case of FIG. 10, the potential is unknown at the node ⊚ and the normal differential of potential is unknown at the node xe2x97xaf, and there are 35 unknown quantities in the whole region VP.
Next, a method for obtaining an unknown quantity in this region VP and its theoretical ground will be discussed. Hereinafter, a region to be analyzed (referred to as xe2x80x9canalysis spacexe2x80x9d) is represented by V and a close boundary existing inside the analysis space V is represented by S. When the analysis space V is assumed to be an electrical element, the inside of the boundary S is conductive and the outside is dielectric.
From the Gauss"" law,
div {right arrow over (D)}=xcfx81xe2x80x83xe2x80x83(1)
is true. Further, from the definition of electric field and dielectric flux density,
{right arrow over (D)}=∈{right arrow over (E)}=xe2x88x92∈ grad "psgr"xe2x80x83xe2x80x83(2)
is also true where {right arrow over (D)}, xcfx81, ∈, {right arrow over (E)} and "psgr" represent a dielectric flux density vector, a charge density, a dielectric constant, an electric field vector and a potential, respectively. An arrow above the reference sign indicates a vector sign.
Substituting Eq. 2 into Eq. 1, the Poisson""s equation as expressed by                               div          ⁡                      (                          grad              ⁢                              xe2x80x83                            ⁢              ψ                        )                          =                  Δψ          =                      -                          ρ              ϵ                                                          (        3        )            
is obtained. If xcfx81=0. Eq.3 is the Laplace""s equation.
The Green""s function G({right arrow over (x)}, {right arrow over (xxe2x80x2)}) is defined as a function satisfying
xcex94G({right arrow over (x)},{right arrow over (xxe2x80x2)})=xe2x88x92xcex4({right arrow over (x)}xe2x88x92{right arrow over (xxe2x80x2)})xe2x80x83xe2x80x83(4)
G({right arrow over (x)},{right arrow over (xxe2x80x2)})=0 (where |xxe2x80x2|xe2x86x92xe2x88x9d)xe2x80x83xe2x80x83(5)
where {right arrow over (x)}, {right arrow over (xxe2x80x2)}, xcex94 and xcex4 represent a position vector for a given position in the analysis space V, a position vector for a given point on the boundary S in the analysis space V, Laplacian with respect to {right arrow over (x)} and Dirac""s delta function, respectively. The Green""s function G({right arrow over (x)}, {right arrow over (xxe2x80x2)}) is a solution to the boundary problem with respect to {right arrow over (x)} when the unit charge is put on the position of {right arrow over (xxe2x80x2)}, from physical interpretation.
The Dirac""s delta function xcex4 satisfies, from its definition,
∫vdVxe2x80x2"psgr"({right arrow over (xxe2x80x2)})xcex4({right arrow over (x)}xe2x88x92{right arrow over (xxe2x80x2)})="psgr"({right arrow over (x)})xe2x80x83xe2x80x83(6)
where dVxe2x80x2 represents a variable of integration in the analysis space V.
Spatial dimension for {right arrow over (x)} and {right arrow over (xxe2x80x2)} is arbitrary. It is known that for example, in a case of two dimensions, the Green""s function G({right arrow over (x)}, {right arrow over (xxe2x80x2)}) is expressed as                               G          ⁡                      (                                          x                →                            ,                                                x                  xe2x80x2                                →                                      )                          =                              -                          1                              2                ⁢                π                                              ⁢          ln          ⁢                      "LeftBracketingBar"                                          x                →                            -                                                x                  xe2x80x2                                →                                      "RightBracketingBar"                                              (        7        )            
and in a case of three dimensions, the Green""s function is expressed as                               G          ⁡                      (                                          x                →                            ,                                                x                  xe2x80x2                                →                                      )                          =                              1                          4              ⁢              π                                ⁢                      1                          "LeftBracketingBar"                                                x                  →                                -                                                      x                    xe2x80x2                                    →                                            "RightBracketingBar"                                                          (        8        )            
Applying the Green""s theorem to the analysis space V and the boundary S using the potential "psgr"({right arrow over (xxe2x80x2)}) on the boundary S and the Green""s function G({right arrow over (x)}, {right arrow over (xxe2x80x2)}),
∫sdSxe2x80x2{right arrow over (nxe2x80x2)}xc2x7{G({right arrow over (x)},{right arrow over (xxe2x80x2)})∇xe2x80x2"psgr"({right arrow over (xxe2x80x2)})xe2x88x92"psgr"({right arrow over (xxe2x80x2)})∇xe2x80x2G({right arrow over (x)},{right arrow over (xxe2x80x2)})}=∫vdVxe2x80x2{G({right arrow over (x)},{right arrow over (xxe2x80x2)})xcex94xe2x80x2"psgr"({right arrow over (xxe2x80x2)})xe2x88x92"psgr"({right arrow over (xxe2x80x2)})xcex94xe2x80x2G({right arrow over (x)},{right arrow over (xxe2x80x2)})}xe2x80x83xe2x80x83(9) 
is true where dSxe2x80x2, {right arrow over (nxe2x80x2)}, ∇xe2x80x2 and xcex94xe2x80x2 represent a variable of integration on the boundary S, an outward unit normal vector on the boundary S, a gradient with respect to {right arrow over (xxe2x80x2)} and Laplacian with respect to {right arrow over (xxe2x80x2)}, respectively.
Transforming Eq. 9 by using Eqs. 3, 4 and 6,
"psgr"({right arrow over (x)})=∫vdVxe2x80x2Gxc2x7({right arrow over (x)},{right arrow over (xxe2x80x2)})xcfx81({right arrow over (xxe2x80x2)})+∫sdSxe2x80x2{right arrow over (nxe2x80x2)}xc2x7{G({right arrow over (x)},{right arrow over (xxe2x80x2)})∇xe2x80x2"psgr"({right arrow over (xxe2x80x2)})xe2x88x92"psgr"({right arrow over (xxe2x80x2)})∇xe2x80x2G({right arrow over (x)},{right arrow over (xxe2x80x2)})}xe2x80x83xe2x80x83(10)
is obtained. In other words, the potential "psgr"({right arrow over (x)}) at a given point (position {right arrow over (x)}) inside the analysis space V can be expressed by using the Green""s function G({right arrow over (x)},{right arrow over (xxe2x80x2)}) and the values of parameters ("psgr"({right arrow over (xxe2x80x2)}), {right arrow over (nxe2x80x2)} and the like) at a point (position {right arrow over (xxe2x80x2)}) on the boundary S.
As to a region having no independent charge therein, such as a dielectric, in the analysis space V, it is assumed that xcfx81=0 in Eq. 10. Then, Eq. 10 is expressed as
"psgr"({right arrow over (x)})=∫sdSxe2x80x2{right arrow over (nxe2x80x2)}xc2x7{G({right arrow over (x)},{right arrow over (xxe2x80x2)})∇xe2x80x2"psgr"({right arrow over (xxe2x80x2)})xe2x88x92"psgr"({right arrow over (xxe2x80x2)})∇xe2x80x2G({right arrow over (x)},{right arrow over (xxe2x80x2)})}xe2x80x83xe2x80x83(11)
In Eq. 11, since the Green""s function G({right arrow over (x)},{right arrow over (xxe2x80x2)}) is known as shown in Eq. 7 or 8 and either one of the potential "psgr"({right arrow over (xxe2x80x2)}) at {right arrow over (xxe2x80x2)} and the normal differential {right arrow over (nxe2x80x2)}xc2x7∇xe2x80x2"psgr"({right arrow over (xxe2x80x2)}) of potential is given as a boundary condition, the other which is unknown can be obtained by performing integral calculus.
When the numerical analysis is performed with computer, approximating the integral, Eq. 11 is expressed in terms of such a matrix equation (where A, xcexa8 and B represent a coefficient matrix, a solution to be obtained and an absolute term, respectively) as
Axcexa8=Bxe2x80x83xe2x80x83(12)
and Eq. 12 is solved numerically to obtain the solution xcexa8.
How to obtain the parasitic capacitance is as follows. From the Gauss"" law, the normal differential {right arrow over (nxe2x80x2)}xc2x7∇xe2x80x2"psgr"({right arrow over (xxe2x80x2)}) of potential is integrated with respect to the variable of integration dSxe2x80x2 on the boundary S between the conductor and the dielectric to obtain the amount of electric charges. With variation of the potential "psgr"j of the j-th conductor, the amount of electric charges Qi of the i-th (ixe2x89xa0j) conductor induced thereby is computed, to obtain the parasitic capacitance Ci,j between the i-th conductor and the j-th conductor as                               C                      i            ,            j                          =                              ∂                          Q              i                                            ∂                          ψ              j                                                          (        13        )            
In these methods, however, an increased number of meshes imposes overload on calculation performance of a CPU and memory capacity in the computer since fragmented nodes and boundary elements are used in accordance with complicated shape of the material to meet the requirement for higher accuracy of analysis, and as a result, it is disadvantageously impossible to broaden an analyzable region. Specifically, the matrix equation of Eq. 12 becomes larger and the use of memory region and calculation time increase. Moreover, it is disadvantageously impossible to obtain approximate solution until the computer completes the matrix equation. These problems especially become prominent in obtaining the potential and parasitic element constant inside a semiconductor device of three-dimensionally complicated shape.
The analysis made by the boundary element method advantageously uses less equations than those made by the finite element method and the finite difference method since the unknown quantity exists only on the boundary. On the other hand, the equations to be solved are disadvantageously complicated. Therefore, the analysis by the boundary element method is not much different from those by the finite element method and the finite difference method in terms of calculation time and the amount of calculations.
In any method, there are still problems that the use of memory region and calculation time increase and the approximate solution can not be obtained until the equations are completely solved. Further, the accuracy of calculated value is found after the equations are completely solved, and therefore it is impossible to perform calculation with desired accuracy determined in advance.
The present invention is directed to a simulation method. According to a first aspect of the present invention, the simulation method comprises the steps of: (a) expressing a first function representing a first physical quantity (potential, temperature) satisfying a Poisson""s equation with respect to a first variable indicating a given position in an analysis space comprising at least one region having a boundary comprising information on position and boundary condition in a form of a first integral with a difference, as a first integrand, between a product of a second function representing the first physical quantity with respect to a second variable indicating a position on the boundary and a first derivative representing a normal component on the boundary out of a differentiated result of a first Green""s function on the first and second variables and the first physical quantity with respect to the second variable and a product of the first Green""s function and a second derivative representing a normal component on the boundary out of a differentiated result of the second function with respect to the second variable, the boundary as a first integral region and the second variable as a first variable of integration; and (b) substituting the first integral sequentially into either one of the second function and the second derivative that is not given as the boundary condition while updating the first variable of integration.
According to a second aspect of the present invention, the simulation method of the first aspect further comprises the steps of: (c) expressing the first function as a sum of series and obtaining values of first to m-th (mxe2x89xa71) terms out of the sum of series by using the information on position and boundary condition that the boundary comprises; and (d) adding all the values of the first to m-th terms obtained in the step (c) to numerically calculate a value of the first function.
According to a third aspect of the present invention, in the simulation method of the first aspect, the at least one region means a plurality of regions, and the method further comprises the steps of: (c) expressing a third function representing a second physical quantity (the amount of electric charges) in a first region which is one of the plurality of regions in a form of a second integral with a product of a third derivative representing a normal component on the boundary of the first region out of a differentiated result of the first function with respect to the first variable on the boundary of the first region and a physical constant depending on a material of the first region as a second integrand, the first variable as a second variable of integration and the boundary of the first region as the second integral region; and (d) performing the step (b) on a fourth function representing a parasitic element constant between the first region and a second region which is one of the plurality of regions obtained by differentiating the third function with respect to the second function in the second region to express the fourth function as a sum of series and obtaining values of first to n-th (nxe2x89xa71) terms out of the sum of series by using the information on position and boundary condition that the boundary comprises; and (e) adding all the values of the first to n-th terms obtained in the step (d) to numerically calculate a value of the fourth function.
According to a fourth aspect of the present invention, the simulation method of the first aspect further comprises the steps of: (c) expressing a third function representing a second physical quantity (current) in the region in a form of a second integral with a product of a third derivative representing a normal component on the boundary out of a differentiated result of a first function with respect to the first variable on the boundary and a first physical constant depending on a material of the region as a second integrand, the first variable as a second variable of integration and the boundary as the second integral region; and (d) performing the step (b) on the third function to express the third function as a first sum of series and obtaining values of first to p-th (pxe2x89xa71) terms out of the first sum of series by using the information on position and boundary condition that the boundary comprises; and (e) adding all the values of the first to p-th terms obtained in the step (d) to numerically calculate a value of the third function.
According to a fifth aspect of the present invention, the simulation method of the fourth aspect further comprises the step of: (f) dividing difference between the second function of a position in the region and the second function of another position in the region by the third function, to numerically calculate a value of a fourth function representing a parasitic element constant (parasitic resistance, parasitic thermal resistance) in the region.
According to a sixth aspect of the present invention, the simulation method of the fourth aspect further comprises the steps of: (f) expressing a fourth function representing a parasitic element constant (self-inductance) in the region in a form of a fourth integral, using a third variable and a fourth variable both of which indicate positions in the region and a fifth function expressed in a form of a third integral with a product obtained by multiplying a fourth derivative representing a differentiated value of the first function with respect to the fourth variable in the region by the first physical constant and a second Green""s function on the third and fourth variables and the first physical quantity as a third integrand, the fourth variable as a third variable of integration and the region as a third integral region, with a product obtained by multiplying an inner product of the fifth function and a fifth derivative representing a differentiated value of the first function with respect to the third variable in the region by the first physical constant, a second physical constant depending on a material of the region and the inverse of the square of a value of the third function as a fourth integrand, the third variable as a fourth variable of integration and the region as a fourth integral region; (g) performing the step (b) on the fourth function to express the fourth function as a second sum of series and obtaining values of first to q-th (qxe2x89xa71) terms out of the second sum of series by using the information on position and boundary condition that the boundary comprises; and (h) adding all the values of the first to q-th terms obtained in the step (g) to numerically calculate a value of the fourth function.
According to a seventh aspect of the present invention, in the simulation method of the first aspect, the at least one region means a plurality of regions, and the method further comprises the steps of: (c) expressing a third function representing a second physical quantity in a first region which is one of the plurality of regions in a form of a second integral with a third variable representing a position on a first boundary which the first region has as a second variable of integration, a product of a third derivative representing a normal component on the first boundary out of a differentiated result of the first function with respect to the third variable on the first boundary and a first physical constant depending on a material of the first region as a second integrand and the first boundary as a second integral region; and (d) performing the step (b) on the third function to express the third function as a first sum of series and obtaining values of first to r-th (rxe2x89xa71) terms out of the first sum of series by using the information on position and boundary condition that the boundary comprises; (e) adding all the values of the first to r-th terms obtained in the step (d) to numerically calculate a value of the third function; (f) expressing a fourth function representing a second physical quantity in a second region which is one of the plurality of regions in a form of a third integral with a fourth variable representing a position on a second boundary which the second region has as a third variable of integration, a product of a fourth derivative representing a normal component on the second boundary out of a differentiated result of the first function with respect to the fourth variable on the second boundary and a second physical constant depending on a material of the second region as a third integrand and the second boundary as a third integral region; and (g) performing the step (b) on the fourth function to express the fourth function as a second sum of series and obtaining values of first to s-th (sxe2x89xa71) terms out of the second sum of series by using the information on position and boundary condition that the boundary comprises; (h) adding all the values of the first to s-th terms obtained in the step (g) to numerically calculate a value of the fourth function; (i) expressing a fifth function representing a parasitic element constant between the first and second regions in a form of a fifth integral, using a fifth variable representing a position in the first region, a sixth variable representing a position in the second region and a sixth function expressed in a form of a fourth integral with a product obtained by multiplying a fifth derivative representing a differentiated result of the first function with respect to the sixth variable in the second region by the second physical constant and a second Green""s function on the fifth and sixth variables and the first physical quantity as a fourth integrand, the sixth variable as a fourth variable of integration and the second region as a fourth integral region, with a product obtained by multiplying an inner product of the sixth function and a sixth derivative representing a differentiated result of the first function with respect to the fifth variable in the first region by the first physical constant, a third physical constant depending on a material of the first region, the inverse of the third function and the inverse of the fourth function as a fifth integrand, the fifth variable as a fifth variable of integration and the first region as a fifth integral region; () performing the step (b) on the fifth function to express the fifth function as a third sum of series and obtaining values of first to t-th (txe2x89xa71) terms out of the third sum of series by using the information on position and boundary condition that the boundary comprises; and (k) adding all the values of the first to t-th terms obtained in the step () to numerically calculate a value of the fifth function.
By using the simulation method of the first aspect of the present invention, either one of the second function and the second derivative that is not given as the boundary condition can be expressed in a form of numerically-calculable function.
By using the simulation method of the second aspect of the present invention, since the first function is obtained as the sum of the values of the first to the m-th terms in an expression where the first function is expanded into the sum of series, the calculation can be performed by the number of the terms in accordance with a predetermined accuracy and stopped before the following term without changing the number of boundary elements for a change of calculation accuracy.
By using the simulation method of the third aspect of the present invention, since the fourth function is obtained as the sum of the values of the first to the n-th terms in an expression where the fourth function is expanded into the sum of series, the calculation can be performed by the number of the terms in accordance with a predetermined accuracy and stopped before the following term without changing the number of boundary elements for a change of calculation accuracy. Further, more accurate value of the fourth function can be obtained in consideration of distribution of physical quantities. Furthermore, the value of the fourth function can be obtained without obtaining the value of the first function.
By using the simulation method of the fourth aspect of the present invention, since the third function is obtained as the sum of the values of the first to the p-th terms in an expression where the third function is expanded into the first sum of series, the calculation can be performed by the number of the terms in accordance with a predetermined accuracy and stopped before the following term without changing the number of boundary elements for a change of calculation accuracy. Further, more accurate value of the third function can be obtained in consideration of distribution of physical quantities. Furthermore, the value of the third function can be obtained without obtaining the value of the first function.
By using the simulation method of the fifth aspect of the present invention, since the fourth function is obtained after the third function is obtained as the sum of the values of the first to the p-th terms in an expression where the third function is expanded into the first sum of series, the calculation can be performed by the number of the terms in accordance with a predetermined accuracy and stopped before the following term without changing the number of boundary elements for a change of calculation accuracy. Further, more accurate value of the fourth function can be obtained in consideration of distribution of physical quantities. Furthermore, the value of the fourth function can be obtained without obtaining the value of the first function.
By using the simulation method of the sixth aspect of the present invention, since the fourth function is obtained as the sum of the values of the first to the q-th terms in an expression where the fourth function is expanded into the second sum of series, the calculation can be performed by the number of the terms in accordance with a predetermined accuracy and stopped before the following term without changing the number of boundary elements for a change of calculation accuracy. Further, more accurate value of the fourth function can be obtained in consideration of distribution of physical quantities. Furthermore, the value of the fourth function can be obtained without obtaining the value of the first function.
By using the simulation method of the seventh aspect of the present invention, since the fifth function is obtained as the sum of the values of the first to the t-th terms in an expression where the fifth function is expanded into the third sum of series, the calculation can be performed by the number of the terms in accordance with a predetermined accuracy and stopped before the following term without changing the number of boundary elements for a change of calculation accuracy. Further, more accurate value of the fifth function can be obtained in consideration of distribution of physical quantities. Furthermore, the value of the fifth function can be obtained without obtaining the value of the first function.
An object of the present invention is to provide a simulation method which allows reduction in use of memory region and calculation time when it is required to quickly obtain an approximate solution even with low calculation accuracy and allows a calculated value of high accuracy to be obtained when high calculation accuracy is required even though the use of memory region and calculation time increase.